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Abstract 

Synchronization of oscillations is a phenomenon prevalent in natural, social, and engineering 
systems. Controlling synchronization of oscillating systems is motivated by a wide range of applications 
from neurological treatment of Parkinson's disease to the design of neurocomputers. In this article, we 
study the control of an ensemble of uncoupled neuron oscillators described by phase models. We examine 
controllability of such a neuron ensemble for various phase models and, furthermore, study the related 
optimal control problems. In particular, by employing Pontryagin's maximum principle, we analytically 
derive optimal controls for spiking single- and two-neuron systems, and analyze the applicability of the 
latter to an ensemble system. Finally, we present a robust computational method for optimal control of 
spiking neurons based on pseudospectral approximations. The methodology developed here is universal 
to the control of general nonlinear phase oscillators. 

Index Terms 

Spiking neurons; Controllability; Optimal control; Lie algebra; Pseudospectral methods. 

I. Introduction 

Natural and engineered systems that consist of ensembles of isolated or interacting nonlinear 
dynamical components have reached levels of complexity that are beyond human comprehen- 
sion. These complex systems often require an optimal hierarchical organization and dynamical 
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Structure, such as synchrony, for normal operation. The synchronization of oscillating systems 
is an important and extensively studied phenomenon in science and engineering [[T]|. Examples 
include neural circuitry in the brain [2J, sleep cycles and metabolic chemical reaction systems in 
biology [|3]|, (H, Q, |l6l, semiconductor lasers in physics [7|, and vibrating systems in mechanical 
engineering Such systems, moreover, are often tremendously large in scale, which poses 
serious theoretical and computational challenges to model, guide, control, or optimize them. 
Developing optimal external waveforms or forcing signals that steer complex systems to desired 
dynamical conditions is of fundamental and practical importance BUl, [fTOll . For example, in neu- 
roscience devising minimum-power external stimuli that synchronize or desynchronize a network 
of coupled or uncoupled neurons is imperative for wide-ranging applications from neurological 
treatment of Parkinson's disease and epilepsy [fTTIl . [fT2ll . lfT3l to design of neurocomputers [14J, 
[fTSl : in biology and chemistry application of optimal waveforms for the entrainment of weakly 
forced oscillators that maximize the locking range or alternatively minimize power for a given 
frequency entrainment range |l9l, [fT6ll is paramount to the time- scale adjustment of the circadian 
system to light [fTTll and of the cardiac system to a pacemaker [18]. 

Mathematical tools are required for describing the complex dynamics of oscillating systems 
in a manner that is both tractable and flexible in design. A promising approach to constructing 
simplified yet accurate models that capture essential overall system properties is through the 
use of phase model reduction, in which an oscillating system with a stable periodic orbit 
is modeled by an equation in a single variable that represents the phase of oscillation [17], 
[fT9l . Phase models have been very effectively used in theoretical, numerical, and more recently 
experimental studies to analyze the collective behavior of networks of oscillators ll20l . [|2T]|. 
[|22l . [|23l . Various phase model-based control theoretic techniques have been proposed to design 
external inputs that drive oscillators to behave in a desired way or to form certain synchronization 
patterns. These include multi-linear feedback control methods for controlling individual phase 
relations between coupled oscillators ll24l and phase model-based feedback approaches for 
efficient control of synchronization patterns in oscillator assemblies ifTOl . [|25l . Il26ll . These 
synchronization engineering methods, though effective, do not explicitly address optimality in 
the control design process. More recently, minimum-power periodic controls that entrain an 
oscillator with an arbitrary phase response curve (PRC) to a desired forcing frequency have 
been derived using techniques from calculus of variations [16]. In this work, furthermore, an 
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efficient computational procedure was developed for optimal control synthesis employing Fourier 
series and Chebyshev polynomials. Minimum-power stimuli with limited amplitude that elicit 
spikes of a single neuron oscillator at specified times have also been analytically calculated 
using Pontryagin's maximum principle, where possible neuron spiking range with respect to the 
bound of the control amplitude has been completely characterized [|27ll . [|28ll . In addition, charge- 
balanced minimum-power controls for spiking a single neuron has been thoroughly studied [29], 

m. 

In this paper, we generalize our previous work on optimal control of a single neuron ll27l . [|28l . 
[|29l to consider the control and synchronization of a collection of neuron oscillators. In particular, 
we investigate the fundamental properties and develop optimal controls for the synchronization 
of such type of large-scale neuron systems. In Section |Ill we briefly introduce the phase model 
for oscillating systems and investigate controllability of an ensemble of uncoupled neurons for 
various phase models characterized by different baseline dynamics and phase response functions. 
Then, in Section |llll we formulate optimal control of spiking neurons as steering problems and in 
particular derive minimum-power and time-optimal controls for single- and two-neuron systems. 
Furthermore, we implement a multidimensional pseudospectral method to find optimal controls 
for spiking an ensemble of neurons which reinforce and augment our analytic results. 

II. Control of Neuron Oscillators 

A. Phase Models 

The dynamics of an oscillator are often described by a set of ordinary differential equations 
that has a stable periodic orbit. Consider a time-invariant system 

X = F{x,u), x{0) = xq, (1) 

where x{t) E M" is the state and u{t) G M is the control, which has an unforced stable attractive 
periodic orbit 7(t) = 7(t + T) homeomorphic to a circle, satisfying 7 = F(7, 0), on the periodic 
orbit T = {y E MJ^ly = 7(t), for < t < T} C M". This system of equations can be reduced 
to a single first order differential equation, which remains valid while the state of the full system 
stays in a neighborhood of its unforced periodic orbit ||3TI . This reduction allows us to represent 
the dynamics of a weakly forced oscillator by a single phase variable that defines the evolution 
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of the oscillation, 



de 

'di 



f{9) + Z{9)u{t), 



(2) 



where 6 is the phase variable, / and Z are real-valued functions, and m G W C M is the external 
stimulus (control) [[3TI . [|32l . The function / represents the system's baseline dynamics and Z 
is known as the phase response curve (PRC), which describes the infinitesimal sensitivity of 
the phase to an external control input. One complete oscillation of the system corresponds to 
9 E [0,27r). In the case of neural oscillators, u represents an external current stimulus and / is 
referred to as the instantaneous oscillation frequency in the absence of any external input, i.e., 
M = 0. As a convention, a neuron is said to spike or fire at time T following a spike at time if 
6{t) evolves from ^^(0) = to 9(T) = 2n, i.e., spikes occur at 6* = 2mT, where n = 0, 1, 2, . . .. 
In the absence of any input u{t) the neuron spikes periodically at its natural frequency, while by 
an appropriate choice of u{t) the spiking time can be advanced or delayed in a desired manner. 

In this article, we study various phase models characterized by different / and Z functions. In 
particular, we investigate the neural inputs that elicit desired spikes for an ensemble of isolated 
neurons with different natural dynamics, e.g., different oscillation frequencies. Fundamental 
questions on the controllability of these neuron systems and the design of optimal inputs that 
spike them arise naturally and will be discussed. 

B. Controllability of Neuron Ensembles 

In this section, we analyze controllability properties of finite collections of neuron oscillators. 
We first consider the Theta neuron model (Type I neurons) which describes both superthreshold 
and subthreshold dynamics near a SNIPER (saddle-node bifurcation of a fixed point on a periodic 
orbit) bifurcation II33H. Il34l. 



1) Theta Neuron Model: The Theta neuron model is characterized by the neuron baseline 
dynamics, f{e) = (1 + /) + (1 - /) cos 9, and the PRC, Z{9) = 1 - cos 61, namely. 



where / is the neuron baseline current. If / > 0, then f{9) > for all t > 0. Therefore, in the 
absence of the input the neuron fires periodically since the free evolution of this neuron system. 




(3) 
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Fig. 1. The Free Evolution of a Tlieta Neuron. Tlie baseline current / — 100, and hence it spikes periodically with angular 
frequency a; = 20 and period To — tt/IO. 



i.e., 9 = f{9), has a periodic orbit 

0{t) 



2tan-i I • ' for / > 0, 



(4) 



with the period Tq = iv/y/l and hence the frequency u = 2vT, where c is a constant depending 
on the initial condition. For example, if 9(0) = 0, then c = 0. Fig. \T\ shows the free evolution 
of a Theta neuron with / = 100. This neuron spikes periodically at Tq = vr/lO with angular 
frequency tu = 2n/T = 20. When / < 0, then the model is excitable, namely, spikes can occur 
with an appropriate input u(t). However, no spikes occur without any input u{t) as 



for / < 0, 



and there are two fixed points (one of which is stable) for u{t) = 0. 

Now we consider spiking a finite collection of neurons with distinct natural oscillation fre- 
quencies and with positive baseline currents. This gives rise to a steering problem of the finite- 
dimensional single-input nonlinear control system, 6 = /(6) + Z{Q)u{t), where G ^2 C M", 
f,Z:Q^ M", and u E U C M. In the vector form, this system appears as 



" 9i ' 




ai + f3i cos 9i 




' 1 


— cos 9i 






a2 + (^2 cos 92 


+ 


1 


— cos 6*2 






an + f3n COS 9n 




1 


— cos 9n 



u{t), 



(5) 
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in which = 1+/^ = l+uf/A, Pi = 1-1^ = l-ujf/A, and > for all i G A/" = {1, 2, . . . , n}. 
Note that ctj + /3j = 2 for all i E J\f. The ultimate proof of our understanding of neural systems 
is reflected in our ability to control them, hence a complete investigation of the controllability 
of oscillator populations is of fundamental importance. We now analyze controllability for the 
system as in ([5]), which determines whether spiking or synchronization of an oscillator ensemble 
by the use of an external stimulus is possible. 

Because the free evolution of each neuron system 9i, i E J\f, in ([5]) is periodic as shown in dU, 
the drift term / causes no difficulty in analyzing controllability. The following theorem provides 
essential machinery for controllability analysis. 

Theorem 1: Consider the nonlinear control system 

x{t) = f{x{t)) + uit)g{x{t)), x{0)=xo. (6) 

Suppose that / and g are vector fields on a manifold M. Suppose that {/, g} meet either of the 
conditions of Chow's theorem, and suppose that for each initial condition xq the solution of 

is periodic with a least period T(xo) < P E Then the reachable set from xq for Q is 
{exp{/, gjLAja^o^ where {/, gjiA denotes the Lie algebra generated by the vector fields / and 
g, and {explf, g}LA}G is the smallest subgroup of the diffeomorphism group, diff(M), which 
contains exptrj for all rj E {/,(?} la ll35l . 

Proof. See Appendix |Al □ 
The underlying idea of this theorem for dealing with the drift is the utilization of periodic 
motions along the drift vector field, f{x), to produce negative drift by forward evolutions for 
long enough time. More details about Theorem [T] can be found in Appendix |Al Having this 
result, we are now able to investigate controllability of a neuron oscillator assembly. 
Theorem 2: Consider the finite-dimensional single-input nonlinear control system 

Q = f{Q) + Z{Q)u{t), (7) 
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where = (^i, ^2, • • • , On)' G ^] C M" and the vector fields /, Z : ^] ^ R" are defined by 





ai + (3i cos 9i 




1 


— COS Q\ 


/(e) = 


02 + 1^2 cos 02 


, z(e) = 


1 


— cos Q2 








Oln + /3n COS 6'„ 




1 


— COS Qn 



in which ai = 1 + li, /Si = 1 — Li, and /j > for all z G A/" = {1, 2, . . . , n}. The system as in 
([7]) is controllable. 

Proof. It is sufficient to consider the case where ^ Ij for i ^ j and i, j G A/", since otherwise 
they present the same neuron system. Because /j > for all i G M, the free evolution, i.e., 
u(t) = 0, of each 9i is periodic for every initial condition ^^(O) G M, as shown in dH), with the 
angular frequency Ui = l^fTi and the period Tj = Tr/v^- Therefore, the free evolution of 6 is 
periodic with a least period or is recurrent (see Remark [I]). We may then apply Theorem [T] in 
computing the reachable set of this system. Let 

adgh{Q) = [g,h]{e) 

denote the Lie bracket of the vector fields g and h, both defined on an open subset of M". 
Then, the recursive operation is denoted as 

adj/i(e) = [g,ad''g-^h]{Q) 

for any k > 1, setting ad^h{Q) = h(Q). The Lie brackets of / and Z include 

(_l)fc-i2*--(ai -/3i)^-isin^i " 
(-l)'=-i2^(a2-/32)'-'sin^^2 

(-l)'=-i2^(ai - f3if~\ai cos 9^ + (5i) ' 
{-lY-^2\a2 - P2Y-\a2 cos ^2 + P2) 

for k G Z+, positive integers. Thus, {/, adJZ}, m G Z+, spans at all G ^2 since a^ — (3i = 
ujf/2 and coi are distinct for i G M. That is, every point in M" can be reached from any initial 



adf-'Z 



adfZ = 
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condition 6(0) G f], hence the system ^ is controllable. Note that if 6 = 0, {adfZ}, k E Z+, 



Remark 1: If there exist integers rrij and rij such that the periods of neuron oscillators are 
related by mjTj = rijTj for all pairs, i,j G N, then the free evolution of is periodic 
with a least period. If, however, such a rational number relation does not hold between any two 
periods, e.g., Ti = 1 and T2 = \/2, it is easy to see that the free evolution of is almost- 
periodic [36] because the free evolution of each 6i, i G A/", is periodic. Hence, the recurrence 
of / in (|7]) together with the Lie algebra rank condition (LARC) described above guarantee the 
controllability. 

Controllability properties for other commonly-used phase models used to describe the dynam- 
ics of neuron or other, e.g., chemical, oscillators can be shown in the same fashion. 

2) SNIPER PRC: The SNIPER phase model is characterized by / = the neuron's natural 
oscillation frequency, and the SNIPER PRC, Z = z{l — cos 9), where z is a model-dependent 
constant [|33l . In the absence of any external input, the neuron spikes periodically with the period 
T = 2n/uj. The SNIPER PRC is derived for neuron models near a SNIPER bifurcation which is 
found for Type I neurons [iMl like the Hindmarsh-Rose model [l37l. Note that the SNIPER PRC 
can be viewed as a special case of the Theta neuron PRC for the baseline current / > 0. This can 
be seen through a bijective coordinate transformation 6'(0) = 2 tan^^lv^ tan(0/2 — tt/2)] + ir, 
(j) G [0, 27r), applied to ©, which yields ^ = w + ^(1 - cos (j))u{t), i.e., the SNIPER PRC with 
z = 2/ CO. The spiking property, namely, 9{(l) = 0) = and ^(0 = 2n) = 2n is preserved under 
the transformation and so is the controllability as analyzed in Section III-B1[ 

More specifically, consider a finite collection of SNIPER neurons with /(0) = (cui, a;2, • • • , i^n)' 
and Z(Q) = {zi(l — cos 9i), ^2(1 — cos 6*2), ... , z„(l — cos 6'„))', where conventionally, Zi = 2/uji 
for i G A/". Similar Lie bracket computations as in the proof of Theorem |2] result in, for 
k = 1, . . . ,n. 



and thus span{/, Z}la = IR", since Ui 7^ Uj for i 7^ j. Therefore, the system of a network of 
SNIPER neurons is controllable. 



spans M". 



□ 



adf-^Z= {{-l)''-^ziujf^-^sm9 
adfZ= cos^i, 



.,i-l)'-'z^ul'-'sm9^)\ 

{-lf-hnUl''cOs9n)', 
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3) Sinusoidal PRC: In this case, we consider /(6) = (wi, U2, ■ ■ ■ , oJn)' and Z{Q) = {zi sin 9i, 
Z2 sin 02, ■ . ■ , Zn sin On)', where > and zi = 2/cOi for i = 1, . . . ,n. This type of PRC's with 
both positive and negative regions can be obtained by periodic orbits near the super critical 
Hopf bifurcation[31 ]. This type of bifurcation occurs for Type II neuron models like Fitzhugh- 
Nagumo model [381 • Controllability of a network of Sinusoidal neurons can be shown by the 
same construction, from which 

adf-'Z = z.iol^'^ cos 01, i-lf-^znool''-' cos^„)', 

adfZ = {{-1)^ZiUjI^ sin ^1, ... , sin^„)', 

and then span{/, Z}la = I^" for coi ^ coj, i ^ j. Therefore the system is controllable. 

III. Optimal Control of Spiking Neurons 

The controllability addressed above guarantees the existence of an input that drives an en- 
semble of oscillators between any desired phase configurations. Practical applications demand 
minimum-power or time-optimal controls that form certain synchronization patterns for a pop- 
ulation of oscillators, which gives rise to an optimal steering problem, 

min J = ip{T,eiT))+ I C{Q{t),u{t))dt 

Jo 

s.t. e{t) = f{e) + z{e)u{t) (8) 
9(0) = Go, e(r) = 

\u{t)\<M, VtG[0,T], 

where G M", m G M; : K x M" M, denoting the terminal cost, £ : M" x M ^ M, denoting 
the running cost, and /, Z : — )■ M" are Lipschitz continuous (over the respective domains) 
with respect to their arguments. For spiking a neuronal population, for example, the goal is to 
drive the system from the initial state, ©q = 0, to a final state 6^ = (2mi7r, 2m27r, . . . , 2m„7r)', 
where rrii G Z+, i = 1, . . . ,n. Steering problems of this kind have been well studied, for example, 
in the context of nonholonomic motion planning and sub-Riemannian geodesic problems [39], 
[|40]| . This class of optimal control problems in principle can be approached by the maximum 
principle, however, in most cases they are analytically intractable especially when the system is 
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of high dimension, e.g., greater than three, and when the control is bounded, i.e., M < oo. In 
the following, we present analytical optimal controls for single- and two-neuron systems and, 
furthermore, develop a robust computational method for solving challenging optimal control 
problems of steering a neuron ensemble. Our numerical method is based on pseudo spectral 
approximations which can be easily extended to consider any topologies of neural networks, e.g., 
arbitrary frequency distributions and coupling strengths between neurons, with various types of 
cost functional. 

A. Minimum-Power Control of a Single Neuron Oscillator 

Designing minimum-power stimuli to elicit spikes of neuron oscillators is of clinical impor- 
tance, such as deep brain stimulation, used for a variety of neurological disorders including 
Parkinson's disease, essential tremor, and Dystonia, and neurological implants of cardiac pace- 
makers, where mild stimulations and low energy consumption are required [fT2l . PTI . Optimal 
controls for spiking a single neuron oscillator can be derived using the maximum principle. In 
order to illustrate the idea, we consider spiking a Theta neuron, described in ([3]), with minimum 
power. In this case, the cost functional is J = u^{t)dt, and the initial and target states are 
and 27r, respectively. We first examine the case when the control is unbounded. 

The control Hamiltonian of this optimal control problem is defined by if = M^ + A(a+/9 cos 6+ 
u — ucosO), where A is the Lagrange multiplier. The necessary conditions for optimality yield 
A = -H = A(/3 - M)sin6', and u = -|A(1 - cos6') by |^ = 0. With these conditions, the 
optimal control problem is then transformed to a boundary value problem, which characterizes 
the optimal trajectories of 9{t) and A(t). We then can derive the optimal feedback law for spiking 
a Theta neuron at the specified time T by solving the resulting boundary value problem. 



More details about the derivations can be found in Appendix IB-1[ 

Now consider the case when the control amplitude is limited, namely, \u(t)\ < M, Vt G [0, T]. 
If the unbounded minimum-power control as in Q satisfies < M for all t E [0, T], then 



-(a + 13 cos 6) + J(a + /3 cos - 2Ao(l - cos 

u {0} = 

1 — cos u 

where Aq = A(0), which can be obtained according to 



(9) 




(10) 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (t/T) 

Time (t/T) 

(a) (b) 

Fig. 2. |(a)| Minimum-power controls, M^(i;), for spildng a Theta neuron withi lj = 1 at various spiking times that are smaller 
(T = 3,4) and greater (T — 8) than the natural spiking time To = 27r subject to M = 1. |(b)| The resulting optimal phase 
trajectories following u^it). 



the amplitude constraint is inactive and obviously the optimal control is given by ^ and (flOl) . 
However, if |m*(6')| > M for some 9 E [0,2n], then the optimal control u%.j is characterized by 
switching between u* (9) and the bound M (see Appendix IB-2|) . 

-M, u*{9) < -M 
u*m{9) = \ u*{9), -M <u*{9) < M (11) 
M, u{9y > M, 

where the parameter Aq for u* [9) in dH) is calculated according to the desired spiking time T by 

Jq a + p cos 9 + [1 — cos 9)u\j 
The detailed derivation of the control u\,i is given in Appendix IB-21 Fig. [2] illustrates the 
optimal controls and the corresponding trajectories for spiking a Theta neuron with natural 
oscillation frequency u; = 1, i.e., / = 0.25, a = 1.25, and (3 = 0.75, at various spiking times 
that are smaller (T = 3,4) and greater (T = 8) than the natural spiking time To = 2tt with the 
control amplitude bound M = 1. Because the unconstrained minimum-power controls for the 
cases T = 4 < To and T = 8 > To, calculated according to satisfy \u*{9)\ < M, there are 
no switchings in these two cases. 
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B. Time-Optimal Control of Two Neuron Oscillators 

Spiking a neuron in minimum time, subject to a given control amplitude, can be solved in a 
straightforward manner. Consider the phase model of a single neuron as in ([21), it is easy to see 
that for a given control bound M > 0, the minimum spiking time is achieved by the bang-bang 
control 

' M, z{e) > 
-M, z{e) < 0, 

which keeps the phase velocity, 9, at its maximum. The minimum spiking time with respect to 
the control bound M, denoted by T*f„, is then given by 



(13) 



M 



1 



eeA 



f{e) + z{e)M 



de + 



fie) - z{e)M 



(14) 



where the sets A and B are defined as 

A = {e\ z{e) > 0,0 < ^ < 27r}, 
B = {e\ z{e) <o,o <e <2n}. 

Time-optimal control of spiking two neurons is more involved, which can be formulated as in 
([8]) with the cost functional J = Idt and with 

e{t) = f{Q) + z{e)u{t), (15) 

where 

fi ai + f3icos9i Zi 1- cos 6^1 

/ = = , Z = = . (16) 

/2 a2 + P2 cos 62 Z2 1 - cos 62 

Our objective is to drive the two-neuron system from the initial state 60 = (0, 0)' to the desired 
final state 6^ = (2mi7r, 2m2vr)' with minimum time, where mi,m2 E Z+. The Hamiltonian for 
this optimal control problem is given by 



H = \o + {\,f + Zu), 



(17) 



where Aq G M and A G are the multipliers that correspond to the Lagrangian and the system 
dynamics, respectively, and ( , ) denotes a scalar product in the Euclidean space E^. 

Proposition 1: The minimum-time control that spikes two Theta neurons simultaneously is 
bang-bang. 



January 5, 2012 



DRAFT 



13 



Proof. The Hamiltonian in (fTTl ) is minimized by the control, 

= { (18) 

1^ -M for (pit) > 0, 

where is the switching function defined by = (A, Z). If there exists no non-zero time interval 
over which = 0, then the optimal control is given by the bang-bang form as in (fTSi) . where 
the control switchings are defined at = 0. We show by contradiction that maintaining = 
is not possible for any non-zero time interval. Suppose that 0(t) = for some non-zero time 
interval, t E [ti, T2] , then we have 

0=(A,Z) = O, (19) 
0=(A, [f,Z]) = 0, (20) 

where [/, Z] denotes the Lie bracket of the vector fields / and Z. According to (fT9l) and (|20l) . 
A is perpendicular to both vectors Z and [f,Z], where 

2 sin 01 
2 sin 6*2 

Since A 7^ by the non-triviality condition of the maximum principle, Z and [/, Z] are linearly 
dependent on t G ['ri,r2]. One can easily show that these two vectors are linearly dependent 
either when 9i = 2nn and G ^1 ^ and 62 = 2nn, or 9i = 62 + 2nn and 62 G M, where 
n E Z. These three families of lines represent the possible paths in the state-space where can 
be vanished for some non-trivial time-interval. Now we show that these are not feasible phase 
trajectories that can be generated by a control. Suppose that {9i{t), 6*2 (t)) = {2mT, a) for some 
r > and for some n eTL, where a G M. We then have Q\{t') = 2 7^ 0, irrespective of any 
control input. Hence, the system is immediately deviated from the line Q\ = 2nn. The same 
reasoning can be used for showing the case of 62 = 2mx. 

Similarly, if (^^i(t), 02{t)) = (a + 2n7r, a) for some r > and for some n G Z, in order for 
the system to remain on the line {Oi{t), 62{t)) = {62 + 2n'K, 62), it requires that 6i{t) = 92{t) for 
t > T. However, this occurs only when 9i = 2mn and 62 = 2{n + m)7r, where m G Z, since 
61—62 = (/i — /2 ) ( 1 — cos 61 ) . Furthermore, staying on these points is impossible with any control 
inputs since for 6i(t) = 2mn and 62{r) = 2(n+m)TT, the phase velocities are 6i(t) = 62{r) = 2, 
which immediately forces the system to be away from these points. Therefore, the system cannot 
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be driven along the path {62 + 2mT, 62) ■ This analysis concludes that = and = do not 
hold simultaneously over a non-trivial time interval. □ 

Now, we construct the bang-bang structure for time-optimal control of this two-neuron system 
and, without loss of generality, let Aq = 1. 

Definition 1: We denote the vector fields corresponding to the constant bang controls u{t) = 
—M and u{t) = M hy X = f — MZ and Y = f + MZ, respectively, and call the respective 
trajectories corresponding to them as X- and Y- trajectories. A concatenation of an X-trajectory 
followed by a F-trajectory is denoted by XY , while the concatenation in the reverse order is 
denoted by YX. 

Due to the bang-bang nature of the time-optimal control for this system, it is sufficient for 
us to calculate the time between consecutive switches, and then the first switching time can be 
determined by the end point constraint. The inter- switching time can be calculated following the 
procedure described in 021, El, [gl. 

Let p and q be consecutive switching points, and let pq be a F-trajectory. Without loss of 
generality, we assume that this trajectory passes through p at time and is at q at time r. Since p 
and q are switching points, the corresponding multipliers vanish against the control vector field 
Z at those points, i.e., 

(A(0),Z(p)) = (A(r),Z(g))=0. (21) 

Assuming that the coordinate of j9 = (6'i, O2)' , our goal is to calculate the switching time, r, in 
terms of 61 and 62- In order to achieve this, we need to compute what the relation (A(r), Z{q)) = 
implies at time 0. This can be obtained by moving the vector Z(q) along the F-trajectory 
backward from q to p through the pushforward of the solution uj(t) of the variational equation 
along the F-trajectory with the terminal condition uj{t) = Z{q) at time r. We denote by e*^(p) 
the value of the F-trajectory at time t that starts at the point p at time and by (e^*^)* the 
backward evolution under the variational equation. Then we have 

a;(0) = (e--^)*u;(r) = (e"-^). = (e^^^), Z(e-^(p)) = (e"-^), Z e-^(p). 

Since the "adjoint equation" of the maximum principle is precisely the adjoint equation to 
the variational equation, it follows that the function t 1— )■ {\{t),uj{t)) is constant along the F- 
trajectory. Therefore, (A(r),Z(g)) = also implies that 

(A(0),a;(0)) = (A(0), (e"-^), Z e-^(p)) = 0. (22) 
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Since A(0) 7^ 0, we know from dlB and ^ that the two vectors Z{p) and {e''^^)^ Z e'^^ {p) 
are linearly dependent. It follows that 



-tY\ rr tY 



Ze^^p), 



(23) 



where 7 is a constant. We make use of a well-known Campbell-Baker-Hausdorff formula [i45l 
to expand {e''^^)^ Z e^^ {p), that is, 



n=0 

A straightforward computation of Lie brackets gives 

adyZ = [F, Z] = [f + MZ, Z] = [/, Z] = 2 



sin 9i 
sin 6*2 



adl.Z = [Y,[Y,Z]] = 2if-AZ), 



where A = diag {2(«i - 2 + M), 2(^2 - 2 + M)}, and furthermore 

a4"+iZ = (-1)"2"(A + M/)"[/, Z], 
ad^^+^Z = (-1)"2"+^(A + M/)"(/ - AZ). 

Consequently, we have 
which is further simplified to 



gTady- ^ 



:tl +M — (M — /3i ) cos 0-1 , A/ — /Si — (a-i + Af) cos fili x^, / , sin . / -1 , n t \ 

' 2(ci-l+M) - + l(c.i-l+M) cos(2TVai-l+M)+ ^^^_^'^^^ sin(2rVai-l+J\/) 



3 + M~(M-/32)cos92 , A/~,32-(a2+M) cos flg cos^2T^/ Q:9-l+M ) + 
2{q2-i+m) + 2rr,n-i + Ml cot,(ZTya2 i--tm)+ 



2(c«2-l+M) 



i92_ 



i(2rVa2-l+M) 



This together with (|23l) yields 



(l-cos 6»2) 
=(l-cosei) 



2(ci-l+M) ^ + 2(c\-l+M) cos(2TVai-l+A/)+ ^^^_^'^^^ sin(2rVai-l+Af) 



M-/3i-(Qi+Af) 



"^^^;^^|]r"^ + "-1(;^^r"^ -s(2.v/^^;^)+^|^ sin{2.V7i;3T+M)J . (24) 



This equation characterizes the inter-switching along the F-trajectory, that is, the next switching 
time T can be calculated given the system starting with (^^i, 62) evolving along the F-trajectory. 
Similarly, the inter- switching along the X-trajectory can be calculated by substituting M with 
-M in (El. 
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Note that the solution to (l24l) is not unique, and some of the solutions may not be optimal, 
which can be discarded in a systematic way. The idea is to identify those possible switching 
points calculated from (l24l) with = that also having the appropriate sign for (p. We focus 
on the case where / and Z are linearly independent, since the case for those being linearly 
dependent restricts the state space to be the curve 

(tti + f3i cos^^i)(l — cos 6*2) = (0:2 + /32 cos 6^2) (1 — COS^^i). 

If / and Z are linearly independent, then [/, Z] can be written as [/, Z] = kif + k^Z, where 
^ _ 2sin^i(l -cos^2) -2sin02(l -cos^i) 

[oLi + /?! COS^^i)(l — cos 6*2) — {Ci2 + /32 COs6'2)(l — COs6'i) ' 

^ 2 sin Q\{a\ + /3i cos ^^i) — 2 sin 6'2(a2 + fi^ 00862) 

^ (ai + Pi cos^^i)(l — cos 6*2) — (0:2 + l32 cos 6^2) (1 — cos^^i) 

As a result, we can write (p = (A, [/, Z]) = ki{X, f) + k2{\, Z). Since we know that at switching 

points (f) = (A, Z) = 0, the Hamiltonian, as in (flTI) . H = and the choice of Aq = 1 makes 

(A,/) = —1. Therefore, at these points, we have = —ki, and the type of switching can be 

determined according to the sign of the function ki. If ki > 0, then it is an X to F switch since 

(p < and hence (p changes its sign from positive to negative passing through the switching 

point, which corresponds to switch the control from — M to M as in (fTSi) . Similarly, if ki < 0, 

then it is a F to X switch. Therefore the next switching time will be the minimum non-zero 

solution to the equation (|24|) that satisfy the above given rule. For example, suppose that the 

system is following a F-trajectory starting with a switching point Pi = {0\i02)' ■ The possible 

inter- switching times {Tjj}, j = 1, . . . with 1 < 2 < • • • < Tj„ can then be calculated 

according to (l24l) based on pj. Thus, the next switching point is p.,. = (6'[,6'2)' = e'^'''"^(pj), 

Tj r = min{rj 1, . . . , Tj „}, such that ki{6\, 62) < 0, which corresponds to an F to X switch. 

Now in order to synthesize a time-optimal control, it remains to compute the first switching 

time and switching point, since the consequent switching sequence can be constructed thereafter 

based on the procedure described above. Given an initial state 60 = (0, 0)', the first switching 

time and point pi will be determined according to the target state, e.g., 0^ = (2mi7r, 2m2vr)', 

where mi, 7/12 G Z+, in such a way that the optimal trajectory follows a bang-bang control 

derived based on pi will reach 6;^. Under this construction, we may end up with a finite number 

of feasible trajectories starting with either X- or F-trajectory, which reach the desired terminal 

state. The minimum time trajectory is then selected among them. 
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1.0 1.87 3.0 3.56 4.0 5.0 5.61 h 271 

Time (t) ( 1^=0.3) 

(a) (b) 



Fig. 3. |(a)| Time optimal control for two Theta neuron system with I\ = 0.3 (oi = 1.3, /3i = 0.7) and I2 = 0.9 (q2 = 
1.9, P2 ~ 0.1) to reach (27r,47r) with the control bounded by M = 0.5 and |(b)| corresponding trajectories. The gray and white 
regions represent where k\ is negative and positive, respectively. 



Fig. [3] illustrates an example of driving two Theta neurons time- optimally from (0, 0)' to 
(27r,47r)' with the control bound M = 0.5, where the natural frequencies of the oscillators are 
ui = 1.1 (Ji = 0.3) and uj2 = 1.9 {I2 = 0.9) corresponding to ai = 1.3, /3i = 0.7 and «2 = 1-9, 
(32 = 0.1. In this example, the time-optimal control has two switches at t = 1.87 and t = 3.56 
and the minimum time is 5.61. 

C. Simultaneous Control of Neuron Ensembles 

The complexity of deriving optimal controls for higher dimensional systems, i.e., more than 
two neurons, grows rapidly, and it makes sense to find out how the control of two neurons 
relates to the control of many. One may wonder whether it is possible to use a (optimal) control 
that spikes two neurons to manipulate an ensemble of neurons whose natural frequencies lie 
between those of the two nominal systems. Of course, if trajectories of the neurons with different 
frequencies have no crossings following a common control input, then the control designed for 
any two neurons guarantees to bound trajectories of all the neurons with their frequencies within 
the range of these two nominal neurons, whose trajectories can then be thought of as the envelope 
of these other neuron trajectories. We now show that this is indeed the case. 

Lemma 1: The trajectories of any two Theta neurons with positive baseline currents following 
a common control input have no crossing points. 
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Proof. Consider two Theta neurons modeled by 

6»'i = (1 + /i) + (1-/i)cos^i + (1-cos0i)m, ^i(O) = 0, (25) 

6'2 = (1 + /2) + (1-/2)cos^2 + (1-cos^2)m, ^2(0) = 0, (26) 

with positive baseline currents, h.h > 0, and assume that uJi < uj2, which implies Ji < I2 
since Jj = i = 1,2. In the absence of any control input, namely, m = 0, it is obvious 
that 6*1 (t) < 92(t) for all t > since Ji < I2. Suppose that 6*1 (t) < 92{t) for t e (0, r) and 
these two phase trajectories meet at time r, i.e., 6'i(r) = ^^2(t). Then, we have 9i(t) — 6*2 (r) = 
(Ji — /2)(1 — cos(6'i(r))) < and the equality holds only when the neurons spike at time r, 
i.e., 6i{t) = 62{r) = 2mr, n G Z+. As a result, ^i(r+) < 6'2(r+), because 6i{t) = 62{t) and 
9i{t) < 6*2 (r), and hence there exist no crossings between the two trajectories 9i{t) and 6*2 (t). 
□ 

Note that the same result as Lemma [U holds and can be shown in the same fashion for both 
Sinusoidal and SNIPER phase models, as described in Section ITI-B2I and Ill-B3[ when the model- 
dependent constant zi > 22 if ^1 < which is in general the case. For example, in the SNIPER 
phase model, z conventionally takes the form 2; = 2/a; as presented in Section ITI-B2I 

This critical observation extremely simplifies the design of external stimuli for spiking a 
neuron ensemble with different oscillation frequencies based on the design for two neurons with 
the extremal frequencies over this ensemble. We illustrate this important result by designing 
optimal controls for two Theta and two Sinusoidal neurons employing the Legendre pseudospec- 
tral method, which will be presented in Section |IVl Fig. |4] shows the optimized controls and 
corresponding trajectories for Theta and Sinusoidal neurons with their frequencies u belonging 
to [0.9,1.1] and [1.0,1.1], respectively. The optimal controls are designed based only on the 
extremal frequencies of these two ranges, i.e., 0.9 and 1.1 for the Theta neuron model and 1.0 
and 1.1 for the Sinusoidal model. 

This design principle greatly reduces the complexity of finding controls to spike a large 
number of neurons. Although the optimal control for two neurons is in general not optimal for 
the others, this method produces a good approximate optimal control. In the next section, we will 
introduce a multivariate pseudospectral computational method for constructing optimal spiking 
or synchronization controls. 
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IV. Computational Optimal Control of Spiking Neuron Networks 

As we move to consider the synthesis of controls for neuron ensembles, the analytic methods 
used in the one and two neuron case become impractical to use. As a result, developing 
computational methods to derive inputs for ensembles of neurons is of particular practical interest. 
We solve the optimal control problem in ([8]) using a modified pseudo spectral method. Global 
polynomials provide accurate approximations in such a method which has shown to be effective 
in the optimal ensemble control of quantum mechanical systems [|46l, [|47l|. [|48l . ||49l . Below 
we outline the main concepts of the pseudo spectral method for optimal control problems and 
then show how it can be extended to consider the ensemble case. 

Spectral methods involve the expansion of functions in terms of orthogonal polynomial basis 
functions on the domain [—1, 1] (similar to Fourier series expansion), facilitating high accuracy 
with relatively few terms lISOl . The pseudospectral method is a spectral collocation method 
in which the differential equation describing the state dynamics is enforced at specific nodes. 
Developed to solve partial differential equations, these methods have been recently adopted to 
solve optimal control problems llSTI . [|52l . [|53l . We focus on Legendre pseudospectral methods 
and consider the transformed optimal control problem on the time domain [—1, 1]. 

The fundamental idea of the Legendre pseudospectral method is to approximate the continuous 
state and control functions, 0(t) and u{t), by N^'^ order Lagrange interpolating polynomials, 
/ArG(t) and I^u^t), based on the Legendre-Gauss-Lobatto (LGL) quadrature nodes, which are 
defined by the union of the endpoints, { — 1,1}, and the roots of the derivative of the A^*'^ order 
Legendre polynomial. Note that the non-uniformity in the distribution of the LGL nodes and the 
high density of nodes near the end points are a key characteristic of pseudospectral discretizations 
by which the Runge phenomenon is effectively suppressed ll54ll . The interpolating approximations 
of the state and control functions, 0(t) and u{t) can be expressed as functions of the Lagrange 
polynomials, ikit), [l55l 



N 




(27) 



k=0 



N 



U 



k=0 
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The derivative of /Ar6(t) at the LGL node tj, j = 0,1, . . . , N, is given by 

k=0 k=0 

where Dj^ are elements of the constant (A^ + 1) x (A^ + 1) differentiation matrix IISOll . Finally, 
the integral cost functional in the optimal control problem dS]) can be approximated by the 
Gauss-Lobatto integration rule, and we ultimately convert the optimal control problem into the 
following finite-dimensional constrained minimization problem 



mm — w 



2 ^ -J -J 

j=0 

N _ J, _ 

s.t. ^D,,e, = - /(e,) + w,z(e 

fc=0 



(28) 



e(-i) = o, 

6(1) = e,, Vje{0,l,...,iV}, 

where 0^ = (2mi7r, 2m2ix, . . . , 2m„7r)', rrii E Z+, i = 1, . . . ,n, is the target state and wj are 
the LGL weights given by 

2 1 

Wo 



' iV(iV + l)(L^(t,)P' 
in which L^v is the A^*^ order Legendre polynomial. Solvers for this type of constrained nonlinear 
programs are readily available and straightforward to implement. 

Remark 2 (Extension to an infinite ensemble of neuron systems): The pseudospectral compu- 
tational method can be readily extended to consider an infinite population of neurons, for instance, 
with the frequency distribution over a closed interval, u E [uJa,cOb] C In such a case, the 
parameterized state function can be approximated by a two-dimensional interpolating polynomial, 
namely, Q(t,co) ~ lNxN^Q{t,(^), based on the LGL nodes in the time t and the frequency cu 
domain. Similarly, the dynamics of the state can be expressed as an algebraic constraint and a 
corresponding minimization problem can be formed [|47l . 

A. Optimized Ensemble Controls 

We can now apply the above methodology to synthesize optimal controls for neuron ensembles. 
Since neurons modeled by the SNIPER PRC are special cases of the Theta neuron, here we 
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Fig. 4. The controls (top) and state trajectories (bottom) of Theta (left) and Sinusoidal (right) PRC neuron models (for a — 1, 
P — 0.1, T — 2tt in i29i). The Theta model is optimized for uj G [0.9, 1.1]. The Sinusoidal model is optimized for cj £ [1.0, 1.1]. 
The gray states correspond to uncontrolled state trajectories, and provide a comparison for the synchrony improvement provided 
by the compensating optimized ensemble control. 

consider Theta and Sinusoidal neuron models. The computational method outlined above permits 
a flexible framework to optimize based on a very general cost functional subject to general 
constraints. We illustrate this by selecting an objective of the type, 

J = a||ed-e(T)f + /3 [ u^{t)dt, (29) 

Jo 

which minimizes the terminal error and input energy with a relative scaling given by the constants 
a and /3. In highly complex problems, such as those given by ensemble systems as described 
in Remark [21 this scaling provides a tunable parameter that determines the trade-off between 
performance and input energy. 

Fig. m shows the optimized controls and corresponding trajectories for Theta and Sinusoidal 
neuron models for a = 1, /3 = 0.1, T = 27r, and to belongs to [0.9, 1.1] and [1.0, 1.1] respectively. 
In this optimization, the controls are optimized over the two neuron systems with extremal 
frequencies, whose trajectories form an envelope, bounding the trajectories of other frequencies 
in between, as described in Section IIII-C[ We are able to design compensating controls for the 
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012345 n/2n 3n/2 2n 



t t 

Fig. 5. The controls and amplitude constrained controls, A < 2.5, upper left and right respectively, of a Sinusoidal PRC neuron 
model driving five frequencies, (lji, aj2, 1^3, 1^4, 1^5) ~ (1,2,3,4,5), to the desired targets 0{T) — (27r, 47r, Gtt, Stt, IOtt) when 
T — 27r — 0.5. These controls yield highly similar state trajectories (left, shown for unconstrained control) and spiking sequences 
(right, shown for constrained control), which corresponds to when the state trajectories cross multiples of 2n. Black coloring 
indicates a controlled state trajectory or spike sequence, whereas gray coloring indicates a trajectory or spike sequence without 
control. 

entire frequency band solely by considering these upper and lower bounding frequencies. The 
controlled (black) and uncontrolled (gray) state trajectories clearly show the improvement in 
simultaneous spiking of the ensemble of neurons. While a bound is necessary to provide in 
practice, the inclusion of the minimum energy term in the cost function serves to regularize the 
control against high amplitude values. 

In Fig. [51 we demonstrate the flexibility of the method to drive multiple Sinusoidal neurons 
to desired targets. In particular we seek to simultaneously spike five frequencies with widely 
dispersed frequency values at a time T different from their natural period. In this figure we 
consider the frequencies {ui,U2,co3, 104,105) = (1,2,3,4,5) and design controls to drive these 
systems to (27r, 47r, Gvr, Svr, lOvr), respectively, at a time T = 2n — 0.5. Controls for minimum 
energy (a = 0, /3 = 1) transfer can be designed for both the unconstrained and amplitude 
constrained cases (shown for A < 2.5). In both cases, the state trajectories and spike sequence 
(shown in the lower half of the figure) follow the same general pattern. The spike train shows 
that the controls are able to advance the firing of each neuron so that all five spike simultaneously 
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Fig. 6. The controls and amplitude constrained controls, A < 2, upper left and right respectively, of a Theta PRC neuron 
model driving five frequencies, oj — 0^2, 1^3, i^4, 1^5) = (1,2,3,4,5), to the desired targets S(r) — (27r, 47r, Btt, Stt, IOtt) 
when r = 27r — 0.5. These controls yield highly similar state trajectories (left, shown for unconstrained control) and spiking 
sequences (right, shown for constrained control), which corresponds to when the state trajectories cross multiples of 27r. Black 
coloring indicates a controlled state trajectory or spike sequence, whereas gray coloring indicates a trajectory or spike sequence 
without control. 

at the desired terminal time. Again the gray coloring indicates uncontrolled trajectories or spike 
trains and offers a comparison of improvement in synchrony. 

Similarly, Fig. [6] provides the same presentation as above for the minimum energy transfer for 
Theta neurons of the same frequencies to the same desired targets. In this case the constrained 
control is limited to A < 2. 

V. Conclusion 

In this paper, we considered the control and synchronization of a neuron ensemble described 
by phase models. We showed that this ensemble system is controllable for various commonly- 
used phase models. We also derived minimum-power and time-optimal controls for single and 
two neuron systems. The development of such optimal controls is of practical importance, for 
example, in therapeutic procedures such as deep brain stimulation for Parkinson's disease and 
cardiac pacemakers for heart disease. In addition, we adopted a computational pseudo spectral 
method for constructing optimal controls that spike neuron ensembles, which demonstrated the 
underlying controllability properties of such neuron systems. The methodology resulting from this 
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work can be applied not only to neuron oscillators but also to any oscillating systems that can be 
represented using similar model reduction techniques such as biological, chemical, electrical, and 
mechanical oscillators. A compelling extension of this work is to consider networks of coupled 
oscillators, whose interactions are characterized by a coupling function, H, acting between each 
pair of oscillators. For example, in the well-known Kuramoto's model, the coupling between the 
(i, j)-pair is characterized by the sinusoidal function of the form H{9i,6j) = sm{6i — 6j) [19J. 
The procedure presented in Theorem |2l can be immediately applied to examine controllability 
of interacting oscillators by investigating the recurrence properties of the vector field f + H, 
and the Lie algebra {f + H, Z}la- Similarly, the pseudo spectral method presented in Section 
ITVl and its extension addressed in Remark [2] can be employed to calculate optimal controls for 
spiking or synchronizing networks of coupled neurons with their natural frequencies varying on 
a continuum. 

Appendix A 
Chow's Theorem 

Theorem 3: (Versions of Chow's Theorem) Let {/i(a;), /2(x), . . . , fm{x)} be a collection of 
vector fields such that the collection /2(x), . . . , f mix)} la is 

a) analytic on an analytic manifold M. Then given any point xq E M, there exists a maximal 
submanifold N C M containing xq such that {exp{xi}}G Xq = {exp{xi}LA}G Xq = N. 

b) C°° on a manifold M with dim (span{ fi{x)} la) constant on M. Then given any 
point Xq E M, there exists a maximal submanifold N C M containing xq such that 
{eyip{xi}}GXo = {exp{xi}LA}GXo = N. 

For more details, please see [[35l . 

Appendix B 
Optimal Control of a Single Theta Neuron 
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1 ) Unbounded Minimum-Power Control of a Theta Neuron: The minimum-power control of 
a single Theta neuron is formulated as 



min / u^{t)dt, 
Jo 

s.t. ^ = a + /3cos6' + (1 - cos6')M(t), 
^(0) = 0, e{T) = 2tt. 
We then can form the control Hamiltonian, 

H = u'^ + X{a + l3cose + u-ucos6), (30) 

where A is the Lagrange multiplier. The necessary conditions for optimality from the maximum 
principle yield 

A = -— - = A(/3-M)sine, (31) 
r)fT 

— = 2m + A(l -COS0) = 0. 
ou 

Thus, the optimal control u satisfies 

u = -^X{l-cos9). (32) 

With (|32l) and (|3T1) . this optimal control problem is transformed to a boundary value problem, 
whose solution characterizes the optimal trajectories, 

e = a + /3cose-^{l-cos9f, (33) 

A = A/3 + y(l-cos^)sin^, (34) 

with boundary conditions ^^(0) = and 9(T) = 2tx, while Aq = A(0) and A(T) are unspecified. 

Additionally, since the Hamiltonian is not explicitly dependent on time, the optimal triple 
{X,9,u) satisfies H{X,6,u) = c, VO < t < T, where c is a constant. Together with (|32|) and 
(1301) . this yields 

\2 

A(a + /3cos^) - — (1 -cos^)2 = c. (35) 

Since 9(0) = 0, c = 2Xq, where Aq is undetermined. The optimal multiplier can be found by 

solving the above quadratic equation (l35l) . which gives 

_ 2(q; + /3cos^) ±2^(a + /3cos^)2 -2Ao(l -cos^) 

(1 — cos 6')"' 
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and then, from (|33l ). the optimal phase trajectory follows 

e = T V(a + /3 cos 9y - 2Ao(l - cos 9). (37) 
Integrating (|37l) . we find the spiking time T in terms of the initial condition Aq, 

T= / d9. (38) 

Jo y/{a + /3 cos 9y - 2Ao(l - cos 9) 

Note that we choose the positive sign in (|37]) . which corresponds to forward phase evolution. 
Therefore, given a desired spiking time T of the neuron, the initial value Aq can be found via 
the one-to-one relation in (|38] ). Consequently, the optimal trajectories of 9 and A can be easily 
computed by evolving (|33] ) and (|34l ) forward in time. Plugging (|36l) into (|32l) , we obtain the 
optimal feedback law for spiking a Theta neuron at time T, 

-(g + /3 cos + v/(a + /3 cos g)^ - 2Ao(l - cos 

u{t) = z ^ . (39) 

1 — cos 9 

2) Bounded Minimum-Power Control of a Theta Neuron: Given the bound M on the control 
amplitude, if \u*(t) \ < M for all t E [0, T], then the amplitude constraint is inactive and obviously 
the bounded minimum-power control is given by (|39| ) and (l38l) . If, however, > M for 

some time interval, e.g., t E [ti,t2] C [0,T], which alternatively corresponds to |m*(6')| > M 
for 9(ti) = 9i, 9(t2) = 92, and 9 E [^^1,6^2] C [0,27r], the amplitude constraint is active and the 
optimal control will depend on M. We first consider u*(9) > M for 9 E [9i, 6*2] and observe in 
this case that u{9) = M is the minimizer of the Hamiltonian H as in (l30l) . since H is convex in 
u. The Hamiltonian for this interval is then given by if = + X{a + 13 cos 9 + M — M cos 9). 
Because, by the maximum principle, if is a constant along the optimal trajectory, the Lagrange 
multiplier A is given by, 

~ a + /3cos^ + M-Mcos^' ^ ' 

which satisfies the adjoint equation (|3TI) . Therefore, u{9) = M is optimal for 9 E [9i,92]. The 
value of the constant c = 2Ao can be determined by applying the initial conditions, ^^(0) = 
and A(0) = Aq to (l30l) . Similarly, we can show that u(t) = —M is optimal when u*{9) < —M 
for some 9 E [9^, 94] C [0, 2tt]. Consequently, the constrained optimal control can be synthesized 
according to ([II]) and (fT2l) . 

Note that the number of time intervals that |^i*(6')| > M defines the number of switches in the 
optimal control law. Specifically, if |tt*(6')| > M for n time intervals, then the optimal control 
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will have 2n switches. Here we consider the simplest case, where the optimal control has only 
two switches, which is actually the only case for the Theta neuron model. As a result, suppose 
that u* (9) > M for only one time interval, and then there are two switching angles 9i and 62 at 
which u*{9i) = u*{02) = M. These two conditions, together with (fT2l) . determine the unknown 
parameters 61, 62, and Aq that characterize the bounded optimal control, m^^, as given in (fTTI) 
for the specified spiking time T. Note that the range of feasible spiking times is determined by 
the bound of the control amplitude M. A complete characterization of possible spiking range 
can be found in ^27l . 
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